library(Mediana)

# 1. Data model ----
sample_size <- c(100, 150, 200)

# 1.1 Outcome parameter set
# 1.1.1 Standard set
outcome1.control <- parameters(prop=0.30)
outcome1.test <- parameters(prop=0.50)

# 1.1.2 Optimal set
# outcome2.control <- parameter(prop = 0.30)
# outcome2.test <- parameter(prop = 0.50)

# 1.1.3 Pessimeistic set
# outcome3.control <- parameter(dispersion = 0.5, mean = 13)
# outcome3.test <- parameter(dispersion = 0.5, mean = 7.8)

# 1.2 Build Data model
trial.dm <- DataModel() + 
	OutcomeDist(outcome.dist = "BinomDist") +  # Binary endopoints
	# SampleSize(sample_size) + 
	Sample(id = "Control",
	       outcome.par = parameters(outcome1.control),
	       sample.size = c(25, 50, 100)) + 
	Sample(id = "Test",
	       outcome.par = parameters(outcome1.test),
	       sample.size = c(50, 100, 200))

# 2. Analysis model ----
# 2.1 Statistic test
trial.am <- AnalysisModel() + 
  Test(id = "Control vs Test",
       samples = samples("Control", "Test"),
       method = "PropTest")  # For "Binary endpoint"

# 3. Evaluation model ----
# 3.1 Power analysis
trial.em <- EvaluationModel() + 
  Criterion(id = "Marginal power",
            method = "MarginalPower",
            tests = tests("Control vs Test"),
            labels = c("Control vs Test"),
            par = parameters(alpha = 0.025))

# 4. Clinical Scenario Evaluation ----
sim_parameters <- SimParameters(n.sims = 10,
                                proc.load = "full",
                                seed = 42938001)

trial.rst <- CSE(trial.dm, trial.am, trial.em, sim_parameters)

# 5. Simulation results ----
# 5.1 Console
trial.rst.summary <- summary(trial.rst)

# 5.2 Presentation model
# trial.pm <- PresentationModel() + 
# 	Project(username = "SimTrial",
# 		title = "Simulation Trial")  # TODO
# 		description = "Clinical Trial Simulation") +  # TODO
# 	Section(by = "outcome.parameter") + 
# 	Table(by = "sample.size") + 
# 	CustomLabel(param = "sample.size",
# 		label = paste0("N=", sample_size)) + 
# 	CustomLabel(param = "outcome.parameter",
# 		label = paste0(c("Pessimist", "Standard", "Optimist")))

# GenerateReport(
# 	presentation.model = trial.pm,
# 	cse.results = trial.rst,
# 	report.filename = "TrialSimReport.docx")

# 6. Get the data generated in the simulation
trial.data <- DataStack(data.model = trial.dm,
                        sim.parameters = sim_parameters)

# i th simulation, j th scenario, k sample
# trial.data$data.set[[i]]$data.scenario[[j]]$sample[[k]]

# Extract simulation results data (data.set)
# sim_data <- ExtractDataStack(data.stack = trial.data,
#                             data.scenario = 1, 
#                             # sample.id = "Control",
#                             simulation.run = 5
#                             )$data.set[[1]]$data.scenario[[1]]


ExtractSimData <- function(trial.data, scenario, run) {
  # Run and Scenario scale
  max_run <- length(trial.data$data.set)
  if (run > max_run) run <- max_run
  
  max_scenario <- length(trial.data$data.set[[scenario]]$data.scenario)
  if (scenario > max_scenario) scenario <- max_scenario
  
  # Extract simulation results data (data.set)
  sim_data <- ExtractDataStack(data.stack = trial.data,
                              data.scenario = scenario, 
                              # sample.id = "Control",
                              simulation.run = run
                              )$data.set[[1]]$data.scenario[[1]]
  
  data1 <- sim_data$sample[[1]]$data$outcome
  data2 <- sim_data$sample[[2]]$data$outcome
  id1 <- rep(sim_data$sample[[1]]$id, length(data1))
  id2 <- rep(sim_data$sample[[2]]$id, length(data2))
  rst <- data.frame(data = c(data1, data2), id = c(id1, id2))
  
  return(rst)
}

# 6. Analysis ----

dat <- ExtractSimData(trial.data, 1, 5)

library(sm)
sm.density.compare(dat$data, dat$id)

boxplot(data ~ id, data = dat)

t.test(data ~ id, data = dat)
